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ABSTRACT 


A finite difference technique utilizing an irregular- triangular mesh and 
the Wielandt inverse- iteration method is used to compute noimial mode slosh- 
ing in spheroidal tanks of eccentricities ranging from zero to 0.8 under 
zero and low gravitational conditions for a contact angle of 5 degrees. The 
results are used to calculate, using a finite Fourier series expansion, 
liquid response to sinusoidal, square- wave, and periodic pulse lateral per- 
turbing accelerations. Reduction of either liquid volume or gravity level 
decreases the fxmdamental sloshing frequency. 
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SUMMAET 


The small- amplitude lateral sloshing of an incompressible inviscid liquid 
partially filling a spheroidal tank, as determined by surface and gravita- 
tional forces, is studied for zero- and low-g conditions. The problem is 
formulated in a curvilinear coordinate system parallel to the equilibrium 
free- surface . The first few normal modes of oscillation having one nodal 
diameter and the corresponding eigenfrequencies are calculated numerically 
for a liquid-tank contact angle of 5 degrees. These modes are computed by a 
finite difference method using an Irregular triangular mesh and a Wielandt 
inverse iteration technique. The forced response to lateral perturbing 
accelerations is developed in terns of the calculated modes using a finite 
Fourier analysis. The calculations are carried out for tanks of eccentricity 
0., 0.5, 0.68, and 0.8 and for liquid volumes ranging from I/8 to 7/8 of a 
full tank. The (dimensionless) Bond number (p^ a, and o 

are the liquid density, steady axial acceleration, tank semi- major axis, and 
surface tension) ranges from zero to 100. At the lower values of B^, the 
equilibrium meniscus may leave dry spots at both the bottom and the top of 
the tank. 

General conclusions include: (a) Q?he fundamental sloshing frequency is 

generally an increasing function of Bond number B^ and liquid volume. The 
fundamental sloshing frequency is zero or near zero for B^ = 0- when the 
equilibrium free- surface intersects the tank wall in a single circle and is 
generally positive for B =0 when the tank wall and free- surface intersect 

Oi 

in two circles, (b) Computation of the response to lateral perturbing accel- 
erations can be effected by use of the finite Fourier analysis, but, where 
appropriate, engineering computations may be more easily made using a spring- 
mass analog for normal mode liquid sloshing. Parameters for this analog are 
presented for first mode sloshing. These are adequate for engineering compu- 
tations when the first term is dominant in the Fourier series for the forced 
motion, (c) The irregular- triangular, finite- difference, Wielandt inverse- 
iteration scheme appears adaptable for this type of problem in a wide variety 
of tank shapes when currently available digital computers are fully utilized. 


INTRODUCTION 


In the recent past, considerable attention has been given to the sloshing of 
liquids in containers of the shapes used in modern propulsion systems, such 
as cylindrical containers with spheroidal ends (sometimes inverted) or simple 
spheroidal ones. Interest has focused more recently on sloshing under low-g 
conditions, where surface tension effects can dominate the behavior of the 
liquid. 

There now exists a large body of literature concerning sloshing under condi- 
tions where surface tension is not important. However, comparatively 
little work has been done on investigating low-g sloshing, and much of 

it has been limited to conditions in which the effects of surface tension 

(3 7) 

are slight so that gravitational forces still essentially dominate. ’ 

The object of this study is to investigate linearized low-gravity sloshing in 
spheroidal tanks of zero to moderate eccentricity. Specific objectives in- 
clude: (a) Determination of the normal lateral sloshing modes and frequencies 

as a function of the axial acceleration level, the surface tension and density 
of the liquid, and the liquid volume in the tank, (b) Determination of the 
forced response of the liquid to lateral perturbing accelerations, (c) Compu- 
tation of the lateral forces and moments acting on the tank. 

The linear analysis employed here describes small amplitude oscillations of 
the liquid about its equilibrium shape and neglects the influence of viscosity. 
(Experiments have shown that viscosity can be safely neglected for determining 
the natural frequency of even quite small scale models. ) It is assumed 
that the fluid properties and the. contact angle of the liquid are constant and 
do not vary dynamically. 
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The frame of reference used for this study is fixed to the container as shown 
in Figure 1. The liquid in tanks of a space vehicle will not be acted upon 
directly by a lateral force that acts on the vehicle, but rather acted upon 
Indirectly by the moving walls of the tank. For convenience, however, the 
tank is considered fixed and the liquid as being influenced directly by a 
lateral perturbing acceleration in the tank-fixed coordinate system. 

The mathematical problem describing the sloshing is a linear boundary value 
problem. In the tank-fixed coordinate system it can be assumed that the liquid 
motion is irrotational so that the governing equation is Laplace's equation, 
the equation of continuity for an ideal liquid. The solution of Laplace's 
equation in the domain occupied by the liquid must satisfy a zero nomial velo- 
city boundary condition at the tank walls and dynamical and kinematic conditions 
at the free-surface. The dynamical condition is obtained from the tinsteady 
form of the Bernoulli equation and the kinematic condition relates the velocity 
potential to the motion of the free-surface. The motion of the free surface 
is constrained to preserve the contact angle constant at the tank wall. 

Previous analysis related to sloshing in spheroidal tanks has been limited to 
the simpler high-g (zero surface tension) case. It was found, even for this 
case, that the use of approximate methods is required, in general, for solving 
the governing equations. (9 AO, 11 ) jn the more complex low-g case studied here, 
approximate methods are also required. To solve the boundary value problem, 
an extension of the numerical technique successfiolly applied previously to 
low-g sloshing in hemi spherically bottomed cylindrical tanks is used. (^, 5 ) 

The numerical solution of the sloshing problem is carried out in three stages: 
(a) The equilibrium meniscus shape is computed for each combination of Bond 
nrunber (dimensionless ratio of gravitationaJ. to surface tension forces) and 
liquid voltime. (b) A suitable finite difference mesh is constructed by 
mapping a cross section of the liquid conformally onto the unit circle, by 
constructing a mesh within the circular domain, and by mapping this mesh 
conformally back into the liquid cross section, (c) The eigenfrequencies and 
eigenfunctions of the discretized system of equations are determined. The 
forced response to lateral perturbing accelerations, as well as other engi- 
neering data, is then developed in terms of the solutions. 
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PROBLEM FORMQIATION 


In this section the mathematical formulation for lateral sloshing in a 
spheroidal container is presented. It closely follows that presented pre- 
viously for a hemi spherically bottomed cylindrical container and^ wherever 
possible, similar notation is used to facilitate comparison.^ ^ 5 ) ^ 

new formulation is required primarily because the one used earlier is based 
upon representing the free surface as a single valued function of the radial 
coordinate in circular cylindrical coordinates, a representation that is 
not general enough for a spheroidal container. The possibility of repre- 
senting a free surface double- valued in those coordinates is allowed for 
here by employing a surface polar nonnal coordinate system to express the 
free -surface boundary conditions. 

Surface Polar Ufonnal Coordinates 

Let the unit vectors associated with the coordinate directions in circular 
cylindrical coordinates (r,0,z) be denoted by the triple and let 

the equations of a meridian of the equilibrium free- surface be 

r = R(s) and z = Z(s) , 

where the parameter s is arc length along the curve measured from the 
lowest point on the surface (see Figure l).** Then the position vector of 
a point on the equilibrium free- surface, a surface of revolution about the 
z-axis, can be represented by 

* 

Reference numbers are superscripted in parentheses. 

** 

Variables appearing in this report are dimensionless, unless otherwise 
stated. They are defined in the Nomenclature List, Appendix A ^ where the 
relationships to the physical dimensioned variables are given. 
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The surface polar normal coordinates are defined relative to the equilibrium 
free- surface by taking and e^ to be unit vectors tangent and normal 

to a meridian. Let point from the liquid into the gas and e point in 

the direction of increasing s so that the unit vectors [e .e^.e ] form a 
right-handed triple. This triple is associated with the coordinate directions 
in surface polar normal coordinates (s,0,T|), where T| is the distance of a 
point from the equilibrium free- surface along e^. In terms of circular 
cylindrical coordinates the unit vectors are 

-s " ^s -r e + R k , 

o s r s — n s “r s ^ 

and conversely ^ 2 ^ 

-r " ^s -s ■ k = e^ + R e .* 

r s s s —n — s — s s — n 


As a sphere of radius T) rolls over the equilibrium free- surf ace , the center 
generates a "parallel" surface. These parallel surfaces have no self- inter- 
sections so long as 


T1 < min (r^, R^) ( 3 ) 

where 

h ' and = Cii/Z3]-^ 

are the principal radii of curvature of the equilibrium free- surf ace . Con- 
sequently, a perturbed free- surface can be described in surface polar normal 
coordinates by the relation 


T1 = H(s,e;t) 


(M 


The subscripts on R and Z denote differentiation. Throughout the 
paper such notation is used to denote differentiation of a dependent 
variable with respect to an independent variable, when the context makes 
it clear that such differentiation is intended. 
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provided that departure from the eciuilibrium. surface T) 0 is small enough 
so that 

max H< min (R^,Rg) • (5) 

(The relation (5) is not a significant restriction for the small amplitude 
sloshing problems considered here.) I^t ijf be the angle from e^ to 
(alternately from k to e ) so that Z = sinijt and R = cosiji, then one 
obtains directly that the position vector of a point can be described in 
surface polar normal coordinates by 

r = r(s,e,Tl) = (R - TIZ^) + (Z + TjR^) k ; (6) 

when 'll = H(s,9;t) the point lies on the perturbed surface; when T] = 0 the 
point lies on the equilibrium surface and (6) reduces to (l). 


The Governing Equation 

The liquid is assumed to be incompressible and inviscid and its motion ir- 
rotational so that the velocity potential satisfies Laplace’s equation 


$ 


rr 



$ 

r 




+ $ = 0 
zz 


(7) 


within the liquid, subject to the boundary condition on the container wall 
that the normal velocity be zero. 


|i = 0 (8) 

an 

on w. The boundary conditions on the free- surface present the major 
complexity and are derived in the next parts of this section. 


Free- Surface Boundary Conditions 

Bernoulli Equation . One of the free- surface boundary conditions is obtained 
from the nonsteady state Bernoulli equation, which in surface polar normal 
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i 


i 

I 

1 coordinates takes the form 


p + B^(Z+HEg) - B^^(R-HZ^) cos9 + (l+B^) 


+ [l|v$-gxrll^ 


+ (0 . rf - lyi^ llrll^] = f(t) 


(9) 


on T| = H(s^6jt)^ where the sign of the velocity potential $ is chosen so 
that the velocity y is 


V 


7$ = $ e + — + $ k 

— r-r r z— 


( 10 ) 


Equations (9) and (lO) are for fluid motion relative to the container in a 
coordinate system fixed to the container. The choice of the function on the 
right hand side of ( 9 ) is arbitrary and, for convenience, is set equal to 
Po^ the static liquid pressure at a fixed point on the equilibrium meniscus* 
The liquid pressure at the interface is related to the gas pressure through 
the surface-tension and the mean curvature of the free- surface by means of 
the relation 

Pg - P = (11) 


Substitution of this relation and for f(t) into ( 9 ) results in 


1+B 


2J4- B^(Z+iffi^) + B^^(R-HZ^) COS 0 - (l+B^) 


Ijvf - Oxrjl^ 


(0 * if - Mf lyi' 


( 12 ) 


= Pg - Pq on Tj = H. 


The form that the mean curvature M in ( 12 ) takes in surface polar normal 
coordinates can be calculated from the relation^^^ ^ that for a surface defined 
by F = 0 


2j4= - div( gradF /| |gradF ||) 


(13) 
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The vector grad F/I Urad Fll is the unit normal to the surfaces F = constant, 
and ( 13 ) states that the negative of its divergence when evaluated on the 
surface F = 0 is twice the mean curvature of that surface. The desired 
form for J4 is thus obtained by taking F = 11 - H(s,0;t) in ( 13 ) and using 
the expressions for the divergence and gradient in surface polar normal co- 
ordinates. 


The standard formulas expressing the gradient of a function F and the 
divergence of a vector G in orthogonal curvilinear coordinates specialize 
in the case of surface polar normal coordinates to 


and for 


grad F 




(14) 


G = gi f 


52 £e , 


to 


dlv G = 


W3 


[(’^2h3gl)< 


(h2h^g2)0 + 


^ ^1^2® 3^11- 


( 15 ) 


where h^, and are the magnitudes of the derivatives of the position 
vector (6) with respect to the coordinates, that is. 


3r 


9r 


Sr 


*^1 -s ~ bs ’ ^2 -0 ■ d0 " ^3 -n 91] 


or 


\ 1 - ^(^ss^s - ^s^ss) = 1 - VRi . 

h2 = R - 11 Zs = RCl-VRg)^ and 


( 16 ) 


Observe that h^ > 0 and hg > 0 imply that Ij < min (R^,Rg), which is 
equivalent to ( 3 ). (The identity 


(Z R - Z R = Z + R 
ss s s ss ss 


ss 
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which follows from the fact that the parameter s is arc length, is useful 
in evaluating ) 

If the perturbed surface is 

F(s,e,T 1 ;t) = T] - H(s,0;t) = 0 , (17) 


then (l4) becomes 


^ " h^ -s ■ -0 -n ^ 


so that 


G = 


grad F 
grad F 


^ “A ^ •'A 

Q -s ■ ® q -n 


where 


1/2 


(18) 


Q = h^hg ll_gr^ f11 = 


Using (15), (17) ^ and (18) to evaluate (13) yields the desired expression 
for J4 , 



on 11 = H(s,0;t) with h^ and h^ given by (16). 


Contact Angle Condition . The end condition for the free- surface boundary 
condition (l2) and (19) is given by specifying the contact angle between 
the free- surface and the container wall. The condition that the contact 
angle remain fixed at the value © during the motion is that 

cos® = ng. * (20) 


along the curve of intersection, where is the unit normal to the free- 

surface pointing from the liquid into the gas and n is the unit inward 

normal to the container wall. The expression for n^. is given by (18) when 
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it is evaluated on the free surface T] = H(s,6jt) at the intersection with 
the wall. If the equations of a meridian of the container wall are 

r = X(t) and z = Y(t) , 


where T is arclength measured along the curve from the lower pole of the 
container, then n is given by the expression 


n 

— w 


= - Y e + X k 
T-r T— 


evaluated at the intersection with the free- surf ace . With respect to the 
unit vectors of surface polar normal coordinates the expression becomes, 
by use of (2), 


n =(-YR +XZ)e +(YZ +XR)e . 

-W '■ T S T s'' -s T s T s'' -n 

Substituting into (20) then yields the desired contact angle condition. 


h^H h h 

. - -p (-Y^R^ + X^Z^) . -i- ^ X^R^) 


( 21 ) 


along the curve(s) of intersection 


R-T1Z=X , Z + 'I1R=Y. 
' s s 


Kinematic Condition . The remaining free- surface boundary condition is the 
kinematic condition, which connects the motion of the free surface with the 
velocity potential. It is derived from the total deri'vative with respect to 
time of the equation defining the free- surf ace , which is 


[ll-H(s,0;t)] = + grad[ Tl-H(s,9;t)] • grad $ = 0 on Tj = H 
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Expressing the gradient in surface polar normal coordinates then yields the 
kinematic condition as 


H = $ - 

"t T1 


§ H 
s s 

h? 




on T] = H 


( 22 ) 


Linearization 

The only nonlinearities of the problem appear in the free- surface boundary 
conditions, (l2), (19)^ (2l), and (22). These can be linearized by assuming 
that the perturbing Bond number is small enough so that H is small 

compared to the principal radii of curvature of the equilibrium free- surface, 
that is. 


“ [z - R -ZE - ¥] (23) 

ss s s ss s 

where (ZR-ZR)=ilf is the curvature of a meridian and Z /R = sini|(/R 
ss s s ss ^s p s' ' 

is the other principal curvature. Neglecting the 0(ir) terms in (l2)(all 
angular velocity terms are of this order) then yields 


2 

i (« Vs " i " [(^sA - " (r) ! 


H 


+ (z R - ZR )+:r^-B^Z-B^, HR +B. R cos9 
' ss s s ss' R a Of s tr . 


- (1 + B^) - X = 0 on T] = 0 , 


where the constant X is equal to p^ ■ twice the mean curvature 

at the point at which p^ is evaluated (ll). The terms 


(Z R - ZR )+~-B Z-X = 0, 
' ss s s ss' R O' ' 


i2h) 


form the equation of the equilibrium free- surface . With these deleted, the 
equation becomes 
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R - Z R + 
ss s s ss 


r 

(W 


H 


( 25 ) 


- B HR + B, R COS0 - (l + B^) § = 0 on T| = 0 

O' s tr Of X 


The end conditions for (24) and (25) are found from (2l). For the equili- 
brium free- surface (24), where T] = 0, (2l) yields the terminal end condition 


cos© = Y Z + X R I 

T S T SJS=S 
T=T.. 


(26a) 


U 

U 


where s and t are the values of s and t, respectively, at the 
u u 

uppermost intersection of the equilibrium free- surface and the container 
wall. 


Z(Su) = Y(t^) , R(s^) = X(t^) . 

If there is sufficient liquid to cover the container bottom, then the inter- 
section at s = s^ and t = is the only one, and the initial end condi- 
tion is the symmetry condition about the axis of revolution, namely 

Z = 0 at s = 0 . ■ (26b) 

s 

If there is insufficient liquid to cover the bottom, then there is another 
(lower) circle of intersection at s = 0 and t = and the initial end 

condition is (26a) at s = 0 and ^ = Tjj as well. (Recall that s is 
zero at the lowest point on the equilibrium free- surface -- at the axis r = 0 
if there is only one circle of intersection; at the lower circle of inter- 
section if there are two. ) 

The terminal end condition for (25) is derived from the expressions for the 
unit normals at the intersection of the container wall and the equilibrium 
free surface, which are 
12 
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n = - sin© e + cos© e 

L— wJ -s — n 

T=T.. 

and 

s=s 

T1=0 

respectively. The unit tangent to the container wall at the intersection 
in the direction of increasing t is also required; it is 




cos© e + sin© e . 

— s — n 


At the intersection of a meridian of the perturbed free surface T] = H(s,6;t) 
with the container wall let the values of s, T], and t be denoted by s^, 

, and t, . Then if x denotes the angle from e to e (alternately 

X JL —1* — T 

from k to n ) so that 

— — 


dn 
— w 

dT 


dT — T 




then the linearized expression for the wall normal is 



j~n - (t,-t ) e 
L— w 1 u dT — tJ 


T=T 


u 


to first order. Similarly, since 


= _ M e 

ds ds — s ^ 

one obtains for the surface normal 

s=s^ s=s s=s 

T)=11^ Tl=TlJ 
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to first order. This latter expression becomes 

UJ = [- ""sSs - r -( V'u) i 

when ( 19 ), with O(H^) terms neglected, is used. 

The relationship between (s^-s^), ® observed from a 

curvilinear triangle obtained when the perturbed surface lies "above" the 
equilibrium free surface, by extending the equilibrium free surface beyond 
the intersection with the container wall (see Figure 2). In the triangle 
with vertices at 

[s^,0] = (x(t^^), Y(t^)) , 

[s^,0] , and 

= (x(t^), yCtj^)) , 

the angle at is 0, the angle at [s^,0] is right, and the lengths 

of the sides are s^ - s , H(st,Q), and t, - t . Hence, to first order 

1 u' ' 1' 1 u 

s, - s = H cot© 

1 u 

and ''’1 " *^u ” H/sin© . 

By substitution into the requirement that 

= [sh • Sw] 

T-Ti 

one then obtains the desired linearized terminal contact angle condition 

lU 
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(27a) 


H sin0 + (-1^ cos© - H = 0 
s 'ds dr 


at s = s, T1 = 0. 

When there is insufficient liquid to cover the tank bottom, then the initial 
contact angle condition required at the lower circle of intersection is 
similarly found to be 


H sin© - (4^ cos© - •^) H = 0 
s 'ds dT"^ 


(27b) 


at s = 0, Tl = 0. 

Finally, the linearized kinematic condition (22) is 

“t ' = i “ ^ ° • ( 28 ) 

The task, then, is to solve the linear boundary value problem (j), (8), (25), 
(27), and (28), where the equilibrium free-surface is defined by (2k), (26)^ 
and the given volume of liquid. In the foregoing, it is assumed that 
^ negative values of are not considered in this study. 

The Eigenvalue (Normal Mode) Problem 

The linear boundary value problem posed above is inhomogeneous because of the 
transverse perturbing acceleration proportional to B in (25). The solu- 
tion can be obtained in two parts. The problem is first made homogeneous by 
setting B^^ to zero and solving the normal mode problem. The result is a 
set of eigenfunctions that may then be used in a Fourier series expansion to 
obtain the response to the transverse perturbing acceleration. 

Normal Modes . Let the periodic time dependence and the angular dependence of 
§ and H for the kth normal mode be 
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( 29 ) 


$ = (r,z) COS0 cos(cUj^t) 

H = (s) COS0 sin(cBj^t) . 

The 0 dependence chosen in (29) corresponds to the modes having one dia- 
metral node excited by the lateral perturbing acceleration cos0. Sub- 
stituting (29) into (25) with = 0 and into (28) gives 


(1 + B„) Vic - " 5 K ) - i “k 

\ s/s R 

H^ = 0 



(z R - Z R )^ + 

fz ) 

s 

f] 

L 

ss s s ss \ 

) 



on T1 = 0 . (30) 


and 


ci). 


j^Hj^(s) = (r,z) = (s,0) = $j^^(s,0) 


n 


n 


The boundary conditions on and are the one on H in ( 27 ) and 


and also. 


= 0 on w : 
k ^ 

n 


$j^(0,z) = 0 and H^(0) = 0 , 


(31) 


if the liquid covers the container bottom. Solution of this problem yields 
a set of eigenfunctions eigenvalues and eigenmodes (Hjj) • 


Response to Perturbing Accelerations . For a sinusoidal perturbation of 
amplitude 


B. = sin(u3 t) , 

tr tr ' o ^ 

the velocity potential of the perturbed motion can be represented as a Fourier 
series 

§ = COS0 cos(aj^t) S Aj^$j^(r,z) , (32) 


16 


LOCKHEED MISSILES 6t SPACE COMPANY 


where the solution to the normal mode problem, lj^(r,z), satisfies (30) and 
( 3 l)> but not (25) and (28). Inserting the series into (25) and (28), the 
free- surface and kinematic boundary conditions for the perturbed motion, and 
combining the results lead to 


(1 * B„) r 

on the equilibrium free-surface, r = R(s) and z = Z(s). Now, R can be 
expressed as an expression in the evaluated on r = R(s) and z = Z(s), 


R = Z [R(s), Z(s)] . ( 33 ) 

k 

It can be shown that the B,nd the corresponding form a biorthogonal 

set with the orthogonality condition. 


S * A I* <1= 

o V 


= 0 
^ 0 


for k ^ m , 
for k = m , 


where s^ is the length of a meridian of the equilibrium free-surface, measured 
from the axis when there is a single circle of intersection and from the lower 
circle when there are two. Consequently, the expression for the Fourier co- 
efficients in (33) take the form 




s 


2 

R ds 


R ds 


S; -A 

Thus, for sinusoidal excitation, the coefficient in (32) is 


( 3 M 






(U 


TT^ 


V 




00 


( 35 ) 


and the velocity potential for the sinusoidally perturbed motion becomes 
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B. COS0 cos(u) t) D. tt) 


1 + B, 


k u), - 

k o 


The procedure for obtaining the solution for any periodic perturbation that 
possesses a Fourier series expansion of the form 


B. = B . E sin(io» t) 
tr tr HI m ' o 


is similar. The solution is 


$ = COS0 E $, (r,z) E A, cos(mo3 t) , 
^ m ^ 


where 


. (m) °k ”tr ° “‘o 

\ - TT^ 2 . „2 2 

k o 


Particular perturbations of interest are the square wave and the periodic 
pulsing accelerations, for which the Fourier coefficients are 




0 , m even 


for the square wave 


, m odd 


f O , m even 

4 / , \ 2 . mdTT , , 

- (-1) sin , m odd 


for the periodic pulse , 


where 6 is the ratio of At, the perturbing pulse width, to T = 
half its period (see Figure 3)* 


The explicit expressions for the normal surface displacement at the vall^ 
the item of principal interest, are then found from (28). They are 


LOCKHEED MISSILES & SPACE COMPANY 


H] 


e=0 ° irrr “2 2 *k ^s=Z(s ) sinusoidal, (38) 

a k=l - 0 )_ n _) u< 


s=s 


u 


r=R(s;;) 


«o„ 


0=0 

s=s 


B. ), ® “ sin(mu) t) 

^ E D„ [ r »,. ], 


1 + > 1 ^ 1 f 2 2 2s-' k -'z=Z(s ) 

O' k=l m=l m(u). - m (O^ ) n u< 


u 


m odd 


r=s(s;;) 


for square wave^ 
(39) 


and 

H]«., 


0=0 

s=s 


u 


Kr k * r ® (- 1 )^ sin ^ Sin(nm) t) 

1 + B ir / 2 2 2x ^ \ ^z=Z(s ) 


(iK)) 


for periodic pulse perturbations. 


Mechanical Analogue . For practical design applications it is convenient to 
express each sloshing mode in terms of the fundamental mode of an equivalent 
spring-mass mechanical analogue. To calculate the parameters of this analogue 
it is necessary to obtain expressions for the lateral force and moment acting 
on the container wall ^ vhich^ in turn^ can be found from the pressure. The 
pressure at the wall is 

^ liquid 

(4l) 

Pg , in the gas , 

where the first expression is the linearized form of the Bernoulli equation 
(9)« When there is sufficient fluid to cover the container bottom, the lateral 
force on the wall is given by the integral 


2tt 




p siroc COS0 XdT d0 


This integral, which is taken over the entire container wall, becomes, upon 
substituting (4l) and discarding terms of higher than first order. 
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2tt t. 


^ S COS0 XdT d0 


2n T +H CSC© 
-. r. u u 


- ^ ^ (X + B^Z) sinx COS0 XdT d0 


2n T 




'o ''o 


The 0 integration yields zero for the first integral and for the part of 
the second between the limits 0 and t^. Thus 


F =-(X+BZ) sinx X esc® \ H cos0 d0 
X a u' u u *^o 

p2TT pT 

-(1 + B)V $^siroc COS0 XdT d0 . 

^ r\ ^ 


For the kth normal mode one obtains, after substituting (29) and the second 
of (30), 

P uj, ^(l+B ) p'^^u -| 

F = Tr[-(X+B^Z^)sinx^^X^csc® + ^ f j^siriXXdTjHj^^sin(cOj^t) . (42) 


k *^o 
nu 


The moment of the lateral force is given by the integral 


2tt Tt 




'o “o 


which becomes simply 


2tt t. 


= - (l - ■^) ^ ^ pXY COS0 siriX dT 


b '^o "^o 


for the spheroidal tank with origin at the tank center (Figure la). 
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X co^ 5 Y sirOC • 


Using the same technique as above, one obtains in general 


2n 


M 


= (l - r (\+B Z ) sirtK X Y csc@ \ H cos6 
' ,2' L' a u' u u j u 

b o 

p2TT pT 


cie 


and for the kth normal mode 


M 




= ^(1-;;|)[(X+B^\) sinx^XJ^ CSC© - ^ ij^sim XYdT] Hj^^sin(u^t) , 


k '"o 

nu 


For the case when there is not sufficient liquid to cover the container "bottom^ 
the same procedure yields the expressions 


r 

\ A ^ Biac^x^csce 




) p''u 

r i *k %u sin(aip) , 


(lt4) 


k ~T I, 

nu Z 


and 


M 


= "(1 - :^)[(4-^B^Z^) 3iri<^X^YjjCsce ^ sinx^X^Y^csce 


(U, ^(l+B ) p*^u 

J— ^ \ slWXYdrJ sin(i«^t) . 

nu '^Z 


(^5) 


The point of action on the container axis of the single force to which 

the mechanical analogue is made equivalent is then the ratio 
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9 


= K /F. 




(i^6) 


where z is measured from the center of the tank. Because the amplitude 

of the force Imposed by a spring-mass oscillator, which is made equal to 

the amplitude of F , is k x, and the amplitude of the energy is 

X, K 

/2, the spring constant of the analogue must be 




(i+7) 


where is the amplitude of the energy of the kth mode of the liquid, 


P u 9$ 

V = •^ (l+B ) \ ^ R ds 

k 2 '' a' J k 9n 


(k8) 


For the mechanical analogue to have the fundamental frequency its mass 

must therefore be 



( 49 ) 
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MUMERICAL ANALYSIS 


The common angular variation of the velocity potential and normal displace- 
ment for the kth normal mode is explicitly exhibited by the factor cos0 in 
( 29 ) and (30)« Consequently, the domain in which the eigenvalue problem 
developed in the preceding section must be solved can be taken as a radial 
section of the liquid within the tank in the plane perpendicular to the 
single diametral node. When the volume of liquid is large enough to insure 
that the bottom is completely covered, the radial section is bounded by the 
axis of revolution, the meridian of the tank, and the meridian of the free- 
surface (Figure la). When the bottom is uncovered, the radial section is a 
lunar region between the meridians of the free- surface and' the tank (Figure 
lb). 

Throughout this section of the report, when no confusion will result, the sub- 
script k attached to H and $ will be suppressed to simplify the notation. 
By using the relation (30)# written as 

u)H(s) = $^(r,z) = $^(s,0) , 

the eigenvalue problem can be reformulated in terms of the kth velocity po- 
tential and its normal derivative on the equilibrium free- surf ace . For the 
bottom covered case the complete restatement follows. The velocity potential 
satisfies 




$ + § ) 

r zz' 



within the liquid subject to the boundary conditions 


1 ( 0 , z) = 0 


(50) 


( 51 ) 
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along the axis within the liquid, 


§ = 0 
n 


on the wetted tank meridian, and 


- 5 (« [*J3)3 * 1 Vs " 3 - [<^ss’>s - ^ *n 

( 53 ) 

= (1+B^) 0)^ $ 

on the equilibrium free- surf ace . The free- surface boundaiy condition (53) 
is subject to the end conditions 

s ' ' 

at the intersection between the free- surface and tank meridians (s = s^^,11 = 0) 


$ = 0 
n 


at the intersection of the free- surface with the axis (s — O). The meridian 
of the equilibrium free- surface is the solution to the differential equation. 


r 2-1 

(Z R - ZR )+^ -BZ-X = 
L' ss s s ss' R J a 


subject to the boundary conditions 


cos© = Y Z + X R_ 1___ 

T S T SJS=S 


Z = 0 
s 


at s = 0, that when rotated about the axis encloses a given volume of liquid 
within the spheroid. 
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When the fluid volume is so small that the bottom is uncovered, the axial 
boundary condition (51) disappears. The end condition ( 55 ) for the free- 
surface boundary condition is replaced by the constant contact angle condition 



sin© - cos® - •^) $ = 0 

'ds dr n 


( 59 ) 


at s = 0, = the lower intersection between the free- surface and tank 

meridians. Finally the constant contact angle condition 


cos© = Y Z + X R 

T S T S. 


Js=0 


T=T, 


(60) 


replaces the symmetry boundary condition (58) for the equilibrium free- surface . 

The computation of the eigenvalues, ci?, the eigenfunctions, $, and the eigen- 
modes, H, falls into three major tasks: 

• Determination of the meniscus shapes defining the domains filled 
by given volumes of liquid. 

• Development of irregular triangular meshes within the radial 
sections of the liquid such that the finite difference approxi- 
mations to the eigenvalue problems are symmetric. 

• Solution of the approximate eigenvalue problems. 

Finally, finite Fourier analyses in terms of the approximate eigensystems 
yield, almost immediately, predictions of the shape of the forced response 
of the liquid to lateral perturbing accelerations. Moreover, the values 
needed to set up an equivalent spring-mass mechanical analog ((44) through (49)) 
are an easy by-product of the computation. Each of the major tasks as well 
as the subsidiary calculations will be discussed in turn. 
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Computation of Meniscus Shapes 

The Initial Value Problems . Combining the differential equation for the 
meridian of the equilibrium free- surface (56) with the relation 


R R 
s ss 


Z Z 
s ss 


= 0 


derived from 



(s is arc length) yields a 2 x 2 system with the solution 

R = - CZ and Z = CR 
ss s ss s 


where 

C=BZ + X- Z/R = dt/ds (61) 

O' s' 

is the curvature of the meridian. Consequently, (56) is equivalent to the 
redundant first order system 

Z = W 
s 

W = CU 
s 

u = - cw (62) 

s 

R = U 
U^+W^ = 1 

in which U and W are the r and z velocity components of a point 
moving with unit velocity along a curve with the function C (6l) as 
curvature. When C is positive, the curve curls to the left as s increases. 

Because only the tangent vectors are involved in the boundary condition (57)# 
the difference between the z-data of the free- surface and the tank may be 
absorbed in the function z = Y(t). Thus the tank may be permitted to slide 
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up and down the axis of symmetry relative to the surface. A convenient 
z-datum for the equilibrium free- surface when the bottom is covered is 

Z(0) = 0 . (63) 

The conditions 

R( 0 ) = 0 , U( 0 ) = 1 , and W(o) = 0 ( 64 ) 

then follow from the boundary condition (58). 

On the solution to the initial value problem (62), (63) and ( 64 ), there is 
at least one point (there may be three) at which the boundary condition (57) 
holds. Let the tank be translated so that the meridian passes through such 
a point. If the segment of the solution up to a point where (57) holds lies 
entirely within the tank, then the two- point boundary value problem (62), 
subject to (57) and (58) with \ to be determined so that the volume has a 
given value V, may be replaced by the problem of choosing \ so that the 
rotation of the segment of the solution of (62), (63), and ( 64 ) between the 
origin and a point where ( 57 ) holds encloses a given volume V. 

To obtain an analogous reformulation in terms of an initial value problem 
when the bottom of the tank is dry, it is convenient to choose the z-datum 
for the tank so that the equation of the tank meridian is 

+ Y^/b^ = 1 . (65) 

Choose an initial point 

R( 0 ) = R^ = X(tjj) and Z(o) = = Y(t^) (66) 

on the lower branch of (65). Re-expressing the contact angle condition (60) 
in the form ij; = x ” ® yields 
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and 


(67) 


u(0) = R (o) = |x cos© + y sin® I 

S L T T J 

w(0) = Z^(0) = ^y^ cos® - sin®] 

as the other initial conditions needed to define a solution of (62) as an 
initial value problem. 

In the family of solutions to the initial value problem (62) starting from 
(66) and (67) and characterized by positive values of X, there is one for 
which the boundary condition (57) holds at the first intersection with the 
ellipse (65). Thus the two- point boundary value problem (62) subject to 
the conditions (57) and (60) with X chosen so that the volume has a given 
value V may be replaced by the problem of choosing X and either R^ or 
so that the volume enclosed by the rotation of the segment of the 
solution of (62), (66), and (67) between (Rj^^Z^) and the point where (57) 
holds encloses a given volume V. 

Qualitative Behavior of the Solutions . In two special cases the solution to 
(62) subject to the initial conditions (63) and (64) can be written down. 

For X = 0, the solution is the straight line Z = 0, for all because 

C(0) = 0 implies C = 0. Clearly the tank can be positioned so that a hori- 
zontal plane makes the required contact angle. 

From (61), it follows that 

R^ = cosi|f/C and Z^ = sin\^/C (68) 

are an alternate form for (56) or (62). For = 0, the circle 

r(i|() = (2 sini|()/X and Z(\j() = 2(l-cosi|()/X (69) 
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with radius 2/1x1 and centered at (r = 0, Z = 2 /X) satisfies ( 68 ). The 
curvature of the meridian is the constant C = X/2. An easy geometric con- 
struction shows that the contact angle condition determines the radius and 
hence X. Thus for = 0, the equilibrium free- surface is a part of a 
sphere so long as the circle ( 69 ) intersects the elliptical cross section 
of the spheroid in two and only two points. For the spherical tank this is 
the 6 ase for any volume; and the planar solution is the limit of the spherical 
solutions as X approaches zero. 

For > 0, it can be shown that C = dijr/ds is a positive and monotone 
increasing function of s for 0 ^ i|t ^ tt. Thus from ( 68 ) it is apparent 
that for X > 0 the solution starting from ijt = 0 (see (64)) lies in the 
first quadrant. Moreover, the R coordinate increases with s to a local 
maximum at i]( = rr/2 and then decreases. The Z coordinate increases with 
s to a local maximum at \|t = tt. Because the curvature C is monotonically 
increasing, ( 62 ) implies that the solution curls more and more tightly to 
the left. In fact beyond i|r = tt, it curls so rapidly that it never crosses 
the axis again but becomes globally a spiral in the first quadrant. See 
Bakker^^^^ for such curves. 

Only in the limiting case, B^ = 0, when the solution is a sphere with con- 
stant curvature, does the solution ever return to the axis again as s 
increases. For B^ > 0, the denominator of the term Z^/r prevents the 
curve from crossing the axis. 

For all sufficiently large X, the solution to ( 62 ) or ( 68 ) starting from 
( 63 ) and (64) is so tightly curled in the interval 0 ^ s tt that 

R(tt) < 2b 

holds. For such cases the tank may slide along the common axis of symmetry 
so that the upper branch of its cross section is tangent to the solution 
near \|t = tt. Lowering the cross section slightly will produce a position in 
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which the contact angle condition holds when ® is small. Decreasing X 
will uncurl the solution, raise the point at which the contact angle condition 
holds, decrease the volume enclosed between the meniscus and the tank. 

Thus for large X, volume is a decreasing function of X. Unfortunately there 
may be, for a fixed X, several points on the solution to (62), (63), and ( 64 ) 
that both satisfy the contact angle condition and the condition that the solu- 
tion lie entirely within the tank. Consequently, volume V is a multivalued 
function of X. However, X is single valued as a function of V. 

Applying the uncurling argument to the two- circle case shows that through 
each initial point there is at least one solution to (62), (66), and 

(67) that satisfies the contact angle condition. What is not clear is whether 
there may be more than one. 

The One Circle Algorithm 

The Basic Procedure . For a fixed Bond number, the meridian of a meniscus that 
has a single circle of intersection with the tank is determined by specifying 
X/2, the curvature at the axis of symmetry when the initial point of the 
meridian is the origin. The shape of a spheroidal tank with equatorial semi- 
axis unity is specified by the eccentricity e of its polar cross section or 
by its polar semi-axis b = (l - e^)^^^. Testing is facilitated by fixing 
the center of the tank at (0,0,b) so that its cross section in the plane of 
the meniscus meridian is the ellipse 

r^ + ((z-b)/b)^ = 1 ( 70 ) 

that passes through the origin and is tangent to the meniscus meridian there. 

At each step in the numerical integration of (62), the basic procedure is to 
slide the tank cross section down along the axis of symmetry, common to the 
ellipse and the meniscus meridian, until the former passes through the point 
just generated on the latter. With the ellipse so positioned there are two 
questions to be resolved. 
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(1) Does the meridian so far developed lie entirely within the 
ellipse? 

(2) Could the contact angle condition he satisfied within the integra- 
tion step just completed? 

If both answers are affirmative, the basic procedure concludes by cubic inverse 
interpolation for the value of s at which the contact angle condition holds 
and by calculating the cosine of the contact angle achieved, the volume V 
included between the adjusted meniscus and the tank, and the derivative 
dV/dX. 

Testing for Inclusion . Were the meridian of a meniscus to intersect the tank 

at its equator and to satisfy the contact angle condition there, 

R = costjt = sin© would hold. So long as 
s 

R > sin© (71) 

s 

holds^ an intersection between the two meridians (tank and meniscus) at which 
the contact angle condition holds must lie on the lower branch of the ellipse. 
When (71) fails, a satisfactory intersection must lie on the upper branch. 

In the latter case, the algorithm computes the distance 

Q = z - Z = b(l + [ 1 -R^]^/^) - Z ' (72) 

along the line r = R between the terminal point of the meniscus (R,Z) and 
the upper branch of the fixed ellipse (70). A negative Q shows that the 
meniscus meridian has escaped from the fixed ellipse and is an indication 
that the integration should be terminated. If the curvature X/2 of the 
meniscus exceeds the curvature of the ellipse at the origin, a nonnegative 
Q suffices to guarantee that the entire meniscus lies within the ellipse. 

This has been the case for all menisci with satisfactory intersections on the 
upper branch within the range of eccentricities. Bond numbers, and volumes 
considered in this report. Consequently, detailed checking for inclusion 
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along the upper branch has not been included in the present version of the 
algorithm. 

For intersections along the lower branch of the ellipse, where (71) holds, 
the situation is more complicated. The algorithm computes the distance 

P = z - Z = b(l - [1-R^]^/^) - z 


along the line r = R between the terminal point of the meniscus and the 
lower branch of the fixed ellipse (70) and maintains a maximum of the values 
of P previously computed, namely 


P 

m 


max 

[0,s-As] 


[P] 


for a positive 
Because P = 0 


As less than the length of the last integration step. 

at s = 0, P ^ 0. 

^ m 


If 


P ^ P 
m 


(73) 


holds, then translating the fixed ellipse downward through a distance P 
will still leave the meniscus point at which P was attained on or below 
the tank meridian. Thus when (73) occurs, the algorithm ceases, testing and 
returns to the meniscus integration. 


If 


P > P 
m 


(7^+) 


holds, then translating the fixed ellipse downward through a distance P will 
put all previously computed meniscus points above the translated tank meridian. 
Provided the mesh spacing is small enough. 


dP/ds = cosi|i(tari)(-tanilf ) 
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will be positive throughout an integration interval that contains a point 
s at which the contact angle condition is satisfied because the inequality 
X > if must hold not only at s but also in a neighborhood thereof. Con- 
sequently, P evaluated at s will also exceed P^. Thus with small mesh 
intervals it suffices to check (jh) at the end of an integration step. 


Testing the Contact Angle Condition . Four functions derived from the relations 

X = if + ® and it = X - ® (75) 


by taking sines and cosines are used to test for the satisfaction of the 
contact angle condition. For example, a point s at which 

Fj^(s) = sirX - (sinif cos0 + cosit sin®) 
or 

F (s) = Y - (P sin® + Z cos®) 

Its s 

has a zero is a point at which the contact angle condition ( 75 ) is satisfied. 
It is an appropriate test function to use when 

dF^/ds = dr/ds 

is large so that the change of sign of F^^ between the ends of an integration 

interval is easily detected numerically. So long as small puddles are not in 

consideration, F^ may be used when 1 > R > holds. There are two 

■i- s 

principles for choice among F^ and its three analogues. First, the deriva- 
tive must be large; second, difficulty in evaluating the function up to the 
final subtraction must be avoided. Each test function derived from the 
relation ® = x - f(f violates one of these two principles. 

Because the sine is monotone increasing in [0,tt/2], the inequality 

= sirX > sin^ = Z^ ( 76 ) 
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is equivalent to the inequality 


X > 

derived from ( 75 ). These inequalities hold not only at the point where the 
contact angle condition is satisfied but also in a full neighborhood of such 
a point. Provided that the integration intervals are so small that a compu- 
ted meniscus point terminating an interval within which the contact angle 
condition is satisfied is within the nei^borhood in which (76) holds, it is 
safe to use the failure of (76) as a criterion for discontinuing the testing 
procedure. (Practical experience shows that this is safe except near the 
equator, for example, when R, Y^, and all exceed . 99 *) The basic pro- 
cedure can be visualized by sliding the tank up and down along the axis of 
symmetry, common to the tank and the meniscus, until the tank meridian passes 
through the current terminal point of the meridian. 

The Strategy . When, for a fixed tank and Bond number, the detail of the 
relation between X and V, the volume enclosed by the meniscus, is still 
not determined, the basic procedure may be used to survey all possible points 
on the meridian at which the contact angle condition holds. For such a 
survey the numerical integration continues until the meridian satisfies one 
of three conditions that guarantee no further points at which the contact 
angle holds : 

(1) The current terminal point lies beyond the line R = 1 . 

(2) The distance Q (72) has become negative. 

(3) The meridian has passed beyond the point where it has a second 
horizontal tangent, iji s n. 

By carrying out such a survey for a suitable set of X, enough points on 
the curve representing X as a function of V can be determined to permit 
iterating toward the X which correspond to a given set of V via secant 
interpolation. When X is so close to the value corresponding to a given 
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V that, in a neighborhood of the goal, V is a single valued function of 
X and |dv/dx| is neither too large nor too small, convergence is speeded 
by choosing the next value of \ via Newtonian extrapolation. 

Very near the equator, X as a function of V resembles an extremely flat 
parabola and large changes in volume correspond to very small changes in X. 
Consequently, determining menisci which end near the equator is numerically 
difficult. Nevertheless, the strategy described above, occasionally refined 
by using double precision arithmetic and extremely fine mesh spacing, with 
some human intervention, successfully determined 195 of the I96 menisci 
attempted in the survey reported here. 

The Two Circle Algorithm 

Here it is convenient to fix the elliptical tank cross section in the position 
of (70). The starting point (66) for the meniscus meridian on the lower branch 
of the ellipse may be specified by giving either = X(Tj^) or = Y(t^). 
Choosing the latter as parameter leads to simpler derivatives to evaluate; 
however, this may not be the best choice. The contact angle prescribes the 
initial direction of the curve via (67). To complete the specification of the 
putative meridian, a guess for the value of X is needed. 

The Initial points on the meridian are computed from the power series expan- 
sion of R and Z as functions of s, X, and Z^ about the' initial point. 
The numerical integration is then straightforward and is continued until the 
vertical distance Q (72) between the solution and the fixed ellipse turns 
negative. The point at which the putative meridian crosses the ellipse is 
determined by inverse cubic interpolation. The cosine of the contact angle 
achieved, the volume enclosed between the meridian and the tank, and the 
derivatives needed for the various alternative improvement schemes are also 
computed. 

The most desirable of the alternates is to improve the guess for Z. and X 

At 

by a two-dimensional Newtonian iteration scheme with the goal of simultaneously 
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satisfying the contact angle condition and the requirement that the volume 
have a prescribed value. Sometimes it worked. More often it diverged, but 
in so doing frequently gave information about the volume and contact angle 
at nearby starting points. Another alternate, which proved useful when the 
contact angle was far from being satisfied, is to fix and to improve X 

until the approximation to the contact angle condition is satisfactory and 
then to try the Newtonian iteration again. 

From the experiments in attempting to develop a systematic way for choosing 
starting values, enough information was gained about the relation between X 
and Z, for each of the thirteen two-circle cases in the survey to permit 
the computation of meniscus shapes of comparable accuracy to those produced 
automatically for the one- circle cases. 

Were there a reason to construct a good two- circle algorithm, it would seem 
that a combination of the second alternate with the secant interpolation used 
for the one- circle case followed by a two-dimensional Newtonian iteration 
when the volume was well bracketed might succeed. 

Variable Mesh Spacing 

In the one-circle case the mesh spacing along the free- surface near the axis 
need not be fine. However, when the surface is hi^ly curved, fine mesh 
spacing is needed near the tank wall to represent the surface accurately. 
Moreover, by making it possible to pack mesh triangles into long thin regions, 
finer mesh near the tank wall has helped produce satisfactory irregular tri- 
angular meshes. 

The penalty term in the discretization error for using unequal mesh spacing 
is proportional to the difference in lengths between adjacent intervals. Thus 
meshes with intervals in a geometric progression minimize the penalty and allow 
a wider mesh spacing at one end of a mesh than at the other without materially 
increasing the number of points or incurring an unusually large discretization 
error at any point. 
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It is a simple exercise to verify that the Lagrangian interpolation coefficients 
on such a mesh can be expressed in terms of the constant ratio between adjacent 
intervals and the length of one of the intervals, say the first. The Adams 
integration methods can be obtained by integrating a Lagrangian interpolation 
polynomial (with equal mesh spacing) over the first interval in the interpola- 
tion for the corrector and over the next interval for the predictor. Conse- 
quently, when Adams type integration formulas are derived from a Lagrangian 
interpolation polynomial on a mesh with intervals in a geometric progression, 
it turns out that the coefficients are constants times the length of the 
leading interval. Thus replacing an equal mesh by one in a geometric pro- 
gression as the basis for an Adams type integration requires the following 
additional work: a multiplication to form the current mesh length and two 

extra multiplications per integration step per equation (assuming a single 
correction), one each for the predictor and corrector sums. 

The computation of the constant coefficients for a geometric Adams method is 
easy, provided that the mesh ratio is near 1, as it should be. Here where 
order four is used, the computation involves solving two three-by- three 
systems accurately. Because this computation is so simple, it is relatively 
easy to change the mesh ratio. To switch from one mesh ratio to another, 
four recomputations of the geometric Adams coefficients are needed. Such a 
procedure facilitates the development of the triangular meshes for two -circle 
cases where it is desirable to pack triangles into the two comers. As to 
accuracy for the cases at hand the geometric Adams methods are only a little 
less satisfactory than the ordinary Adams methods. 

Mesh Generation 

The second major task is the development of a suitable irregular triangular 
mesh within half an axial cross section of the liquid- filled region. The 
goal is to obtain a mesh nearly free of obtuse triangles — none near the 
free- surface or upper tank wall. In a regular triangular mesh, lay off along 
a horizontal line the number of mesh points to be used on the equilibrium 
free- surface; at the left end point draw a line downward to the left; at the 
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right end point leave just one triangle in the corner and again draw a line 
downward to the leftj finally close the figure by another horizontal line. 
This parallelogram is the logical diagram; the mesh lines in the interior 
and on the boundary are mapped into curved lines in the interior and on the 
boundary of the half cross section of the liquid- filled region (logical into 
real free- surf ace ) by a numerical algorithm. In addition to the number of 
mesh points on the free- surface and the number of mesh lines parallel to the 
free- surface, the mesh spacing along the free- surface may be varied to force 
more triangles into the corner near the contact angle. The details of the 
mesh are controlled by prescribing the location of some of the physical 
images of the logical boundary points. The conputing mesh points are the 
intersections of the curved mesh lines within and on the boundaries of the 
half cross section. 

Figure 4 shows an irregular triangular mesh used in this study. The free- 
surface and the tank wall, from the axis to the point where only two triangle 
meet at the wall, correspond to the longer, horizontal sides of the parallelo 
gram; the axis and the remainder of the wall correspond to the slanting sides 
The darker lines in the Interior are the images of an integral, rectangular 
coordinate system used to locate vertices in the regular triangular mesh. By 
mentally distorting the mesh so that all triangles are equilateral, the user 
can recover a picture of the underlying logical diagram. 

The ellipse with eccentricity e and unit major semi-axis 

r^ + (z/hf = 1 (77) 

in the £ = r + iz plane is mapped onto the unit circle in the w = x + iy 
plane with the origin going into the origin and the real axis onto the real 
axis by the complex function 

w = k^'^^ sn j^(2K/Tr) arc sin(C/e)J (78) 
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where k is the modulus and K the real quarter period of the Jacobian 
elliptic function sn. The terminology here follows the lilBS 

Handbook of Mathematical Functions. The inverse map 


C = e 


sin 






( 79 ) 


Which carries the unit circle back into the ellipse involves an elliptic 
integral of the first kind. The evaluations of the elliptic functions and 
integrals needed to effect (78) and ( 79 ) are made by algorithms derived from 


those of Bulirsch 
Wynn.^^°^ 


(19) 


and the computation of the complex inverse sine follows 


Using (78) to map the meridian of the free- surface 


C(s) = R(s) + i Z(s) 

into the unit circle produces, from the half cross section of the liquid 

filled region bounded by a part of an ellipse, a standardized domain bounded 

by a piece of the unit circle. Figure 5 shows the standardized domain for 

Figure A triangular mesh is produced within the standardized domain bv 

(21) 

using an algorithm described by Winslow. ' The computing mesh (Figure 4 ) 
is generated by applying the inverse map (79) to the standardized map 
(Figure 5). 


An obvious advantage of the standardized domains is the ease with which points 
and distances on the unit circle may be described in terms of a central angle. 
The decisive advantage is that experience gained in producing a suitable mesh 
in one standard domain may be immediately applied to the development of meshes 
within similar standardized domains. Thus meshes developed for spherical 
tanks, where no mapping is necessary, have frequently served as time-saving 
models of meshes for spheroidal tanks with moderate eccentricity. Because 
(79) preserves angles locally but not in the larger scale of the mesh tri- 
angles, a mesh that is satisfactory in the standardized domain may not be so 
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when mapped into the original elliptical domain. The converse may occur. 
Consequently^ the practical procedure is to contirol the mesh by prescribing 
the location of some boundary points on the standardized domain and to 
decide the acceptability of the mesh in terms of the number and location of 
the obtuse triangles that appear in the final computing mesh. 

As an example consider Figure 6 , the computing mesh, and Figure 7, the 
standardized domain, for a typical two- circle case. In the lower tip of 
Figure 7 there is a cluster of obtuse triangles near the free- surf ace , 
marked by "O” at their centroids. These disappear in the computing mesh. 

Figure 6. 

The logical domain for a two- circle case is an isosceles trapezoid built 
from equilateral triangles with its longer side corresponding to the free- 
surface. The obtuse triangles that are marked on Figure 6 are associated 
with the two comers of the trapezoid where only two triangles meet at a 
boundary point. In Figure 4 there is exactly one such obtuse triangle in a 
similar location. Further adjusting of the input parameters may remove such 
triangles. However, experience shows that, the expenditure of effort is un- 
warranted. Improvement comes slowly with change and more objectionable obtuse 
triangles at the free- surface frequently appear as a consequence of the readjust- 
ment. 

Figure 6 shows a discernable discontinuity of the mesh spacing along the tank 
boundary at the equator. Following the broken line parallel to the darker 
broken lines from the equatorial point on the tank boundary to the free- 
surface locates the equatorial point there. Close examination of the mesh 
spacing along the free- surface will indicate that it decreases in both 
directions from the equator and in part accounts for the desirable crowding 
of mesh triangles into the tips of the half cross section. The broken line 
connecting the images of the equatorial points in the logical domain de- 
composes the isosceles trapezoid into two unequal "right" trapezoids which 

(4) 

resemble the logical diagrams used in the previous report. The user 
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prepares input for each half of the logical diagram separately. This input 
is rotated by tt/2 to the left in a special input program to allow main 
mesh generation routines to operate unchanged, and a special output program 
rerotates the mesh by tt/2 to the right to the position Illustrated. 

The procedures used to control the mesh within the standardized domains are 

(4) 

described in the previous report. They include prescribing the location 
of certain control (mesh) points on the unit circle, the number of mesh 
intervals in the group to be constructed between control points, and the 
length of the first interval of such a group. Whenever the data permits, 
the length of the mesh intervals in a group form a geometric progression; 
otherwise, they are equally spaced. In addition to the parallelogram used for 
Figure 4 and the isosceles trapezoid used for Figure 6, the program permits 

(4) 

using the right trapezoid as in the previous report. Moreover, such 
domains may be stacked on top of one another, provided that the right edge 
of the stack is along a mesh line slanting downward to the left and provided 
the left edge has no segments along horizontal mesh lines. 

All meshes required in this survey can, in principle, be generated by the 
procedure described above. Additional programming is still needed to give 
more control over the location of mesh points along the unit circle to produce 
satisfactory meshes for some domains that are characterized by having free- 
surface and tank meridians of nearly equal length and relatively small center- 
line depth. One difficulty to be overcome is already exhibited in Figures 4 
and 5* The fine mesh along the upper part of the circular boundary in Figure 
5 maps a more widely spaced mesh along the periphery of the ellipse. 

Solution of the Approximate Eigenvalue Problem 

The Matrix Approximations . The introduction of surface polar normal coordi- 
nates results in a new form of the free-surface boundary condition (53) and 
new and easily interpretable constant contact angle conditions (54) and (59). 
With 
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0 = + (1/E^) - [(Z^^E^ - + (Z^/e)2] 


(53) niay be written as 


- 5 ‘n ■ 


Let t, and t„ be mid- points of two adjacent mesh intervals [s. -i^s.] 

1 d. J 

and [s.^s.^, 3 * Integrating (80) with respect to R ds over the interval 
[ti^tg] yields the "balance" equation 


1 1 *2 ^ 

5 Q( = ) V ds . ( 1 -B„) $ 

Js=t^ Js=t^ i 


$ R ds 


and applying the ‘usual two- point approximations for the derivatives and making 
one-point approximations for the unknown functions in the integrals yield 


the finite difference equation 


R(t_) R(t^ ) f} 


R(t ) / R(t ) R(t ) r ”2 \ 

- $ ^ + $ ( ^ + — + \ Q(s) R ds) 

n.,^ s_T-s. n. \s_T-s. s.-s. ^ J4. / 

j+1 j+1 j j \ j+1 j J J-1 ‘^t, / 


R(t, ) ? r ^ 

. $ i__ = |.( 1 +B ) O) \ 

"J -1 ^j-^J -1 J “ \ 


with the two integrals to be approximated during the integration of the free- 
surface. Applying the same argument over an interval [t^, s^], with = s^, 
adjacent to the tank wall yields in view of (56) the finite difference 
equation 
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as the analogue to (8l) incorporating the constant contact angle boundary 
condition. For the tvo-circle case, applying the same argument over an 
interval with s^ = 0, adjacent to the tank wall yields a finite 

difference equation similar to ( 82 ) but incorporating the constant contact 
angle condition (59) • Observe that all of the quantities needed to evaluate 
the difference equations (8l) and ( 92 ) can be obtained in the course of the 
integration of the free- surface. 

Assembling the difference equations and dividing by (l+B^) yields a 
matrix- vector equation 


XA^ = 0 


( 83 ) 


where 

2 and cp^ are vectors of values approximating the potential and its 
normal derivative at mesh points in the free- surf ace j 

2 

X is an approximation to u) j 

A is a diagonal matrix whose entries are 1/2 tt times the 

areas of the zones used to approximate the free surface; and 


T is a tridiagonal, symmetric, irreducible matrix that approxi- 
mates the differential operator on $ in the free- surface 

n 

boundary condition. 


The matrix T is positive definite for all cases in this study with 

^ it is positive definite for all but one of the two-circle cases 

studied with = 0; and it fails to be positive definite for all but one 

of the one-circle cases studied with B =0. 

a 
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The reminder of the mtrix approximtions are described in detail in the 
previous report. The only difference occurs in the two- circle case; the 
condition $ = 0 on the axis is replaced by the condition 
tank wall at one end of each mesh line. However, a brief sketch of the 
arguments is included here. 

At each interior mesh point the operator (50) has a seven- point finite 
difference approximtion connecting the potential values at six neighboring 
points to the value at the central mesh point; at nonaxial boundary points 
a similar approximtion using only points within the closed domin can be 

(21 22 23) 

derived, provided the norml derivative is known. ^ ^ On the free 

surface, (83) provides the missing norml derivatives. Let D be the 

mtrix approximting (50) with $ = 0 on the axis and $ = 0 on the wall 

-1 -1 ^ 

and the free- surf ace . Set S = 2AT A. Partition the vector of potential 
values as C£ = components of values of ^ not on 

the free- surf ace . When D is partitioned conformlly with the finite 

difference approximtion to the noiml mode problem takes the form 



which emphasizes the role of the free- surface . 

For meshes consisting entirely of acute triangles, the mtrix D is symmetric, 
irreducible, positive on the diagonal, nonpositive off the diagonal, diagonally 
dominant and has at most seven nonzero entries in each row. Violent deviations 
from acuteness change the properties of D and affect the final results mrkedly; 
smll departures from strict acuteness away from the free- surface have little 
effect. The matrix S inherits its properties from T; it is tridiagonal, 
symmetric, irreducible and for most cases in this study positive definite. 

Because (8^^) can be reduced to a generalized eigenvalue problem of the form 
[a - XB] £ = 0, standard theory guarantees, when B is positive definite, 
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that (8i|-) has solutions ^ ~ where n is the number of 

points on the free- surface. 

When S fails to be positive definite, the numerical evidence so far 
obtained and the form of S support the conjecture that it has precisely 
one nonpositive eigenvalue and that the others are positive. If this 
is the case, working in the space orthogonal to the eigenvector corresponding 
to the nonpositive eigenvalue should restore the property of the preceding 
paragraph with n reduced by one. 

The Modified Wlelandt Inverse Iteration . Factor the matrix S into left 

and right triangular factors, S = LR, by a Cholesky decomposition modified 

when the product of diagonal entries is negative to set the left diagonal to 

the negative and the right diagonal to the positive value of the square root 

of the absolute value of the offending product. When S is positive definite, 

T 

no modification is necessary and R = L . When S fails to be positive 
definite here, the modification takes place only for the final diagonal entries 

Multiply (84) on the left by 


I 0 


I 0 ’ 


I 0 

0 R. 

y insert I = 

0 L 


0 L-1. 


between the matrix and the vector in (84), and carry out the matrix multipli- 
cations to get an equivalent linear system 


1 

1 — 1 

CVJ 

p 

CVJ 

CVJ 

p 


’z' 

LrD^2 ^11^ " 


. z 


= 0 


(85 


where £ = ^2 _z = L As in the usual Wielandt inverse Iteration^^^^, 

the one modified to take advantage of the form of (85) starts from a guessed 
eigenvalue and a guessed eigenvector The latter is improved 

by solving the linear system 
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( 86 ) 


P 22 

^21^' 

CM 

1 — 1 

_g 

RD^^L-\I 



The guess for the eigenvalue is improved by computing the Rayleigh quotient 



and and replace and z^°^ as the starting values for 

the next step in the iteration. The form of the linear system (86) accounts 
for the speed of the iteration. Because the leading entries in the right 
side are zero, when the system is solved by a factorization method, only a 
portion of the downsweep need be carried out. Moreover, because all the non- 
zero entries in RD,_ are in the right portion of the matrix product, the up- 
12 . (l) 

sweep may be discontinued as soon as the corresponding values of ^ have 

been computed. 

The iteration continues until the ratio some fixed 1, has 

1 1 4-ir 

grown so large that no further improvement can he obtained (lO is the value 
used here) or until the relative error ^ is so small that 

further improvement is unlikely (lO ^ seems suitable for the present problems). 


When the iteration converges^ the approximate potential throughout the tank 

^ is computed to obtain the values along the tank wall needed to evaluate 
the lateral force and moment integrals. The approximate eigenmode H is 
produced by solving the linear equation 

R A H = (X^/^2) z (88) 

and the approximate potential on the free- surface is determined by the matrix 
multiplication 


= L _z , 
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When S is not positive definite and X is negative, the square root in 
(88) is set to one because the eigenmode makes no sense. 

Problems vith B^, = 0 . On physical grounds the fundamental eigenvalue should 
be zero for = 0 when the tank is spherical. A finite- difference approxi- 
mation to a continuous problem can at best yield an approximation to a zero 
eigenvalue. It can be shown for the one- circle problems that zero is not an 
eigenvalue of the approximate problem. Here the first approximate eigenvalue 
is small and positive for V = l/8 and small and negative for V = 7 / 8 . For 
V = 1/2 the approximation seems to be "too good"; the program quits under 
circumstances that suggest that the lowest eigenvalue of the finite differ- 
ence approximation is zero to within machine accuracy. Unfortunately time 
has not permitted modifications to the program that would verify this hypo- 
thesis. 


By continuity, the fundamental eigenvalue is expected to be small for the 
other spheroidal tanks. For all but one of the single circle cases with 
B^ ~ 0 tried in the nonspherical tanks, the program produced an approximation 
for the fundamental eigenvalue. These values are small and uniformly negative. 
The exceptional case seems to be another example of "too good" an approxima- 
tion, like V = 1/2 in the spherical tank. In two- circle cases with B =0 

a 

the program produced complete eigensystems; in one case, V = l/^, e = 0.5, 
the fundamental eigenvalue turned out to be small and negative and in another, 
V = 5/6, e = 0 . 8 , very small and positive. At a fixed volume, as the 
eccentricity of the tank increases, the fundamental eigenvalue also increases. 

The pressing and unsolved problem here is to tell which of the small funda- 
mental eigenvalues are truly nonzero. As eigenvalues of the finite differ- 
ence problem approximating the continuous problefa, they are only slightly 
less accurate than those in the remainder of the study, which are uniformly 
numerically acceptable. The question, therefore, would seem to center on the 
discretization error committed in passing from the continuous to the finite 
difference problem. 


47 


LOCKHEED MISSILES & SPACE COMPANY 


In principle, the program can and sometimes does produce some of the higher 
eigenvalues in the one- circle cases with = 0. When it fails to do so, 
it consistently slips hack into the fundamental— even when provided with 
extremely large guesses for the eigenvalues. This circumstance suggests 
that modifying the construction of the starting eigenvector used in the 
Wielandt inverse iteration would improve the performance of the program for 
B = 0. The fact that the chief difference between the one- and two- circle 
programs lies in this area reinforces this view. 

Evaluation of the Forced Response 

In principle, if enough frequencies and corresponding eigenmodes Hj^(s) 

are known, the shape assumed by the perturbed free- surface H(s,0jt) in response 
to a periodic lateral perturbing acceleration which has a Fourier series 
expansion of the form 


= B , Z C sin(mo) t) 
tr tr m m o 


can be calculated by evaluating the series 

COS0B, C 


(89) 


0 ), -m u) 
k o 


that follows from (37) with D, given by (3^)* 


For sinusoidal perturbations, where = 1 and = 0 for m > 1, the 

inner sum in ( 89 ) disappears and the evaluation of the sum over k is straight' 

forward provided that 00 ^. The numerical problem is to be sure that 

truncating the sum to the small set of natural frequencies and eigenmodes 

available does not omit a physically significant' part of the response. For 

other perturbations, the evaluation of the inner sum presents an additional 

practical problem. The terms of largest magnitude in the sum over m usually 

occur for ■values of m near (u, /cu , and this ratio grows rapidly with k. 

k' o 

A feasible procedure for this initial survey is to carry each sum so far 
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beyond the neighborhood of that there is no question of the convergence 

of the inner sum. 

Thus the main question about the results here is whether enough terms have 
been included in the outer sum (89). With only 25 mesh intervals on the free 
surface, the accuracy of the representation of higher eigenmodes may be 
questioned. Consequently, five eigenmodes were used for the majority of the 
cases. The choice of five was also based on the fact that, with kk mesh 
intervals on the free surface, doubling the number of eigenmodes from four 
to eight (and halving the time step) revealed no significant change in the 
computed response. 

Granted that H(s,0;t) can be evaluated satisfactorily, the question arises 

of when and where the maximum excursion produced by a given forcing function 

appears. An approximation to the desired answer can be given by computing 

H(s,Ojt) for a large number of values of t and noting the time at which 

|H(s,Qt)| achieves its maximum among the values computed. In the present 

survey 30 time steps within 3/8 of the period of the fundamental or 60 time 

steps within 1/2 of the period of the fundamental are used. The shape of 

the perturbed surface is plotted at the time of the observed maximum and at 

intervals of 5 time steps going away from the maximum in both directions. 

The period of the fundamental was chosen as the unit of time to give a 

uniform time scale for varying forcing frequencies u) . 

o 

Because H(s,0;t) represents the deviation of the perturbed surface from the 
equilibritan free-surface in the normal direction, this deviation is added 
vectorially to the free-surface to produce a graph of the perturbed surface 
shape in the plane 0 = 0 at each time plotted (cf. Figure I9). To save 
the user the trouble of reflecting the shape in the equilibrium free-surface, 
the reflected shape in the plane 0 = n is also plotted on the same graph. 

The plotted values are normalized so that the maximum amplitude of the observed 
response ((89 computed with cos9/(l+B^) set equal to l) is equal to 

1/10 of the length of the free-surface (s^ = s^^). Thus the amplitude of the 
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transverse Bond number reported in the plots is the value computed from the 
relation 




s /lo = 6 
max' 


tr 


max 

s^t 


lH(s,0;t)l . 


This normalization is adopted to give the plots a large enough scale to make 

meaningful detail visible and at the same time prevent normals from crossing 

within the distance s /lO of the free- surface (see (3) and (23)). Thus 

max 

the free- surface can be surrounded by a grid consisting of four parallels at 

distances of . l/lO and jt l/20 of s • In addition the normals through 

max 

every fifth point on the free- surface starting from s ^^^ are drawn and 
labeled. 


The question of whether or not the tank meridian should also be included in 
the plot can be argued either way. VJhen it is present, it clutters the graph 
but does serve to emphasize that the present study is limited to small 
amplitude analysis and that the scale of the plots in the normal direction is 
arbitrarily, set large to reveal detail of the computed perturbed shape. 
Clearly the perturbed surface cannot cross the tank wall and must extend to 
the wall when the displacement is positive at s 

The present program can certainly be used (with caution) to determine whether 
the computation of the perturbed shape is of enough practical use to warrant 
further refinement. 

Evaluation of the Mechanical Analog Parameters 
Because 


\ + B Z = 2f4 
01 

and 

siriX = Y = bX/(l-eV)^/^ , 

T 
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the lateral force for the two-circle case (44) may be rewritten as 


with 


and 


+ U)(l+B^) [(5X^)/(l-eSc^)^/^] dxj H(s^)j 

^je = • [x csc®/(l-e^X^)l 

-'s=o ^ 

csc0/(l-e^^)]^ 

S-Su u 


If Ca is set to zero, (90) represents F for the one-circle case (42) 

^ X 

as well. 


The coefficients and must be computed at the beginning and end of 

the final Integration of the equilibrium free- surf ace . A natural form for the 
approximation to the integral in (90) is 

^ [($X^)/(l-eS^^)^/^] dT =Ej Wj (91) 

"^JL 

where J is the set of indices that pick out values of the potential on the 
tank wall from the vector of potential values at all the mesh points. The 
set J and the corresponding weights w^ must be determined while the entries 
in the matrix D are being assembled. However, the final evaluation of the 
integral must wait until the potential is generated. 

Because the tank wall is approximated by straight lines between mesh points, 
a mid- point approximation seems appropriate. Let (X^,Y^) and (X2,Y2) be two 
adjacent mesh points on the tank wall with and as corresponding 

potential values. Set 
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and 


X = (x^ + X_)/2 
m ' 1 2 ' 

At = [(Xg - \f + (Yg - Y^)^] ^ . 

Then the approximation 

' * ‘v> dr/2(l-e\2)^/" 

contributes the weight 

X ^ AT/2(l-e^ (92) 

m ' m ' 

to each of the weights w, and w in (9l)« Each w , except at the ends 
of the tank meridian^ is the sum of two terms of the form (92)* An analogous 
procedure is used to compute the moment (43) or (45). 

T -1 

Because the potential is normalized so that $ S $ = 1, it turns out that 

V = tt/2 (1+B^) ^ R ds = tt(1+B^) w^/h . 

Thus the evaluation of the remaining mechanical analog parameters (46), (47), 
(48), and (49) is straightforward. 
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RESULTS AMD CONCLUSIONS 


The Method 

The problem of low gravity sloshing in spheroidal tanks has previously 
been investigated by approximate mathematical methods and by experimental 
methods. Even when surface tension forces are negligible, approximate 
solutions of the governing equations are the only ones possible. 

No attempts have been made thus far to apply these methods to the case in 
which surface tension forces dominate the liquid motion. Experimental work 
has, with one exception, been limited to the case of zero surface tension, 
and even there physical limitations prevented experiments under conditions 
where surface tension actually dominated the liquid behavior. This report 
is the first attenpt to study lateral sloshing in spheroidal tanks under 
such conditions. 


The highly curved free- surface shapes characterized by large departures 
from the horizontal plane are accurately computed here by numerical means 
properly accounting for the equilibrium contact angle and, by iteration, 
the liquid volume in the spheroidal tank. Solution of Laplace *s equation 
in the domain occupied by the liquid subject to appropriate boimdary con- 
ditions at the tank wall and the free surface is acconqolished using a 
finite difference technique on an irregular triangular mesh. Such a mesh 
is space filling allowing accurate approximation of curved tank and free- 
surface boundaries. The basic restriction inherent in the finite-difference 
approximation is that the contact angle be nonzero. 


Demonstration that the computer program can feasibly produce the desired 
numerical calculations is a major result of this work. The basic program 
used in this project is an adaptation of one previously developed^ for 
calculating low-g lateral sloshing in hemi spherically bottomed cylindrical 
tanks. The principal output from the program includes: 
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2 

1. Normal mode eigenvalues ; 

2. Normal mode eigenranctions ; 

5. Normal eigenmode s ; and 

4. Fourier coefficients, D, . 

These have been used to compute, as a Fourier series expansion, the response 
to sinusoidal, square wave, and periodic pulse lateral perturbing accelera- 
tions and to compute the lateral force and moment Imposed on the tank by 
the liquid in kth mode lateral sloshing. Where appropriate, the force and 
moment can be used to compute parameters of an eqtiivalent spring-mass system 
to facilitate engineering computations. 

The Present Study 

The data reported here is a survey of small amplitude linearized sloshing 
in spheroidal tanks of eccentricity 0, 0.5^ 0.68, and 0.8 with axial Bond 
mjimber ranging from zero to 100 and relative liquid volume ranging from 

l/8 to 7/8 of the tank volume. All calculations have been carried out for 

(4) 

a fixed contact angle of 5 degrees consistent with previous work, ' which 
allows adequate representation of small-contact-angle low-g liquid sloshing 
behavior while keeping numerical difficulties inherent in the method for 
very small contact angles minimal. 

Figures 8-l4 show the meniscus shapes considered in the study. At the higher 
Bond numbers, the liquid cross section ranges from a flat puddle in the tank 
bottom to a deep body enclosing a flattened bubble at the top of the tank. 

At the lower Bond number, the liquid free-surface extends from the lower 
into the upper hemispheroid; and, at s-ufficiently small liquid volume and 
Bond number, the bottom of the tank is uncovered, the liquid being found in 
an annular region around the equator. 

2 

Eigenvalues u) ^. 

The variation of first mode (fundamental) eigenvalues as a function of Bond 
n\miber B^ and liquid volume is shown in Figure 15 and Table I. The first 

p 

mode eigenvalue generally increases with Bond nmber and liquid volvune 

for a given tank shape. The behavior of ^ observed here, as a function 
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of liquid volume for = 100 and e = 0 is consistant with the known 
behavior of cD]_^ in a spherical tank for ~ ” (the case of zero 
surface tension). For k 5 an<i e = 0 , the values of (x>^ in Table I 
replace the values reported earlier (^5) 

for cases in which the liquid lies 

entirely within the lower hemisphere. The new values are consistently higher 
and are more reliable because they are based upon the exact constant contact 
angle condition ((27a) or (5^)) rather than upon the approximation used in the 
earlier computations for the hemisphere. For B^ ^ 2 and e = 0, the effect 
of the approximation in the hemisphere appears to be negligible. (Because the 
constant contact angle condition used on the vertical wall was exact, the data 
reported earlier (^>5) need be replaced by the present values only for B ^ 5 
and e = 0 , cases with large Bond ntnnber and with small fill levels in which 
the liquid lies entirely within the hemispherical bottom. ) The fundamental 
eigenvalue for B^ = 0 , e = 0 is zero. (That it is exactly zero can 

be deduced from physical considerations. The coiBputational results presented 
in Table I are, subject to the discretization error, considered to be consis- 
tent with this value. ) 

2 

Table I also indicates that od^ is very nearly zero in spheroids with e ^ 0 
when the meniscus intersects the tank wall in only one circle and in some cases, 
when it intersects it in two. In other two-circle cases, characterized by small 
fill level and large eccentricity, is well away from zero. Note also that 

first mode eigenvalues for the two-circle cases do not fall on the curves in 
Figure 15. This may be explained as follows. Physically, lateral slosing in 
the one-circle and two-circle cases is much different. In one-circle cases, the 
motion is essentially lateral from one side to the tank to the other; the eigen- 
modes are all odd. First mode sloshing in two-circle cases is circumferential 
with liquid moving from one side of the tank to the other circumferentially; 
higher modes are characterized by an almost vertical motion in each side of the 
tank, and both odd and even modes (i.e. an odd or even number of equilibrium 
free-surface crossings) must be considered. 

The first five eignevalues for all cases considered are listed in Table 

I. Note that all eigenvalues for a given tank shape - Bond number combina- . 
tion increase with liquid volume V except for the two-circle of intersection 
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cases where the higher eigenvalues decrease as liquid volume increases 


Dimensional sloshing frequencies u)j^ can be obtained from values in Figure 
15 or Table I from the relation 


u),. 


= (1)^ ((l+B^)a/pa3) 


1/2 


The centerline depth of liquid under the meniscus and other meniscus geometry 
may be obtained from Figures 8-l4. 


Eigenmode s 

The normal departure of the free surface from its equilibrium shape is shown 

in Figures 16-19 for six representative cases. Low Bond ntunber and high Bond 

number shapes are shown in Figures 12 and 15 to illustrate the influence of 

tank shape on the eigenmodes. It is observed that the first eigenmode for the 

smaller Bond numbers is generally convex up in the region 0 :S s < s (H-i ) , 

-^max 

where H-, is the maximum value of Ht ; whereas, by contrast, H-. is 
Tnax 

generally convex down in down in this same region for Bond numbers greater 
than 10. The dependence on the first eigenmode shape is not otherwise 
(except for two-circle cases) greatly affected by changes in tank shape, 
liquid volume, or Bond number. 


Response to Lateral Perturbations 


Practical application of the results of the present analysis centers on the 
response of the liquid to lateral perturbing accelerations. The response to 
sinusoidal perturbing accelerations is given by 


H 


®tr 
0=0 " 1+B 


a 


sinu) t 
o 


Z 

k=l 




2 


U) 


\U) 


(93) 


where Dj^, defined by (5^), is the Fourier coefficient in the expansion 

K = S Vk- 

For all cases in the survey, the quantities needed to evaluate (93) at 
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s (at the intersection of the tank wall with the equilibrium free surfhce) 

are tabulated for the first five modes: in Table I, Dj^ in II, and 

%(su) in III. (For cases in which the liquid lies entirely within the 
hemispherical bottom, the present tables should be used, as pointed out above, 
rather than those reported earlier (^^5).) 

The forced response to sinusoidal, square wave, and periodic pulse lateral 
perturbing accelerations has been calculated for representative cases. In each 
case the transverse Bond number has been chosen to yield a maximum value of H 
equal to i/io of = Sjj^g^ . Surface -normal coordinates (s,H) where s is 
arclength along the free-surface and H is the displacement normal to the free- 
surface are used in plotting (see Figure I9). Qn each plot a surface -normal 
coordinate grid surrounds the free-surface. The lines of constant H in the 
grid are the parallels to the free-surface at H = ± s^^^lO and H= ± s^^x/SO 
to show where the response reaches the maximum and half the maximum excursion. 
The lines of constant s are the normals to the free-surface drawn through 
s = 0 , s = smax , and every fifth mesh point starting from s ^ay . The 
normals are labeled with their s-coordinates. 

The computed response of the liquid to a sinusoidal, perturbing acceleration 
for the two-circle case = 1, V = 3/8, and e = 0.8 is shown in Figure 
19* The response shown has a maximum excursion from equilibrium of l/lO of 
^u‘ magnitude of the transverse Bond niunber required to force the 

response to have this amplitude is given. The simple form of (93) allows 

easy computation of the response to sinusoidal lateral perturbing accelera- 
tions of different magnitudes and at different ratios of perturbing to 

fundamental frequency. The plot of eigenmodes for this case is also included. 

The response of the liquid to any periodic lateral perturbing acceleration 
is given by 


B 


H 


tr 


0=0 1+B 


a k=l 




00 C sinmo) t 
V m o 

^ 2 2 2 
m=l (0, -m (0 
k o 


( 94 ) 
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where the C 's are the coefficients of the Fourier series expansion of the 
lateral acceleration. The response of liquid to square wave and periodic 
pulse perturbing accelerations for the same case is shown in Figures 20 and 
21 for a forcing frequency 7/l0 the fundamental eigenfrequency. Four plots 
are presented in each case showing the motion in the first half cycle leading 
to an amplitude approximating the maximum. In the last plot for each case, 
the maximum departure of the free- surface from its equilibrium shape is l/lO 
the length of a meridian of the free- surf ace . The magnitude of the transverse 
Bond number B, resulting in this displacement is noted on each plot. The 
response near the maximum to sinusoidal, square wave and periodic pulse per- 
turbations is given for the case V = 3/^^ e = 0 in 

Figure 22, and similarly for the case = 5^ V = 3/^> and e = 0.8 in 
Figure 23 and for B^ = 2, V = 3/8, and e = 0.68 in Figure 24. 

The convergence of the series expansion for the forced response to sinusoidal 
lateral perturbations depends mainly on the requirement that the forcing 
frequency u)^ be markedly different from the natural frequencies For 

u) < (I), , the growth of with k enforces the convergence and the five 

terms from Tables II and III are generally sufficient to give three signifi- 
cant figure accuracy. 

More limitations must be imposed on the evaluation of (94). The inner series 

involving the C 's is clearly convergent, depending only on the ratio 

u) However, the do m inant term may occur far out in the series if 

cUj^/(u^ closely approximates an integer found in the definition of the 

C 's. (in the case of square wave and periodic pvilse perturbations, only 
m 

odd integers are important in this connection. ) Moreover, resonance can 
occur if is a submultiple of u)^. It may be concluded that this method 

of finding the response to lateral perturbing accelerations must be used 
with caution in each specific case so that difficulties arising from sub- 
raultiple resonances may be identified and avoided. 
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The Mechanical Analog 

Mach emphasis has in the past been placed on finding simple approximations 
to the lateral sloshing of liquids in tanks such as equivalent pendulums or 
spring-mass oscillators which, in principle, iirpose the same force on the 
tank as does the sloshing liquid. Such models are also useful in physically 
interpreting the results of studies like the present work. Mechanical analog 
parameters for first mode sloshing are shown in Figures 25-28. The lateral 
force obtained by integrating the con^jonent of the pressxire force acting 
horizontally in the 6=0 direction is presented in Figure 25 . The spring 
constant of the spring-mass oscillator producing lateral forces equivalent 
to first mode sloshing is shown in Figure 26. The lateral force action point 
is shown in Figure 27. This quantity gives the required attachment point 
of the spring-mass oscillator to produce the proper moment about the center 
of the tank. 

Note in Figure 25 that the lateral force imposed by the liquid as it moves 
to the right (i.e. in the 6=0 direction) is negative for smaller values 
This result arises from the fact that the liquid pressTire near the 
free-sxu'face is less than that over the liquid. The integral of this pressure 
depression is embodied in the first term on the right hand side of (^2). It 
is observed that it is negative (i.e. it acts in the 6 = tt direction) and 
offsets the inertial behavior of the liquid embodied in the second term in 
(^2). As the Bond number increases, the first term diminishes, compared to 
the second. When = <», only the second term found in similar analyses 
restricted to the case of zero surface tension remains. For the higher 

values of Bond number considered in this study the lateral force is always 
positive. 

The spring constant and the mass are always positive as a consequence of 
( 47 ) and ( 49 ). Because of this, the spring constant is plotted only in the 
region where the lateral force is positive. The mechanical analog is not 
considered valid when a positive displacement of the right hand 11 mh of the 

free- surface results in a negative force (i.e. one directed to the left in 
Figure l). 
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The magnitude of the mass in the mechanical analog may be obtained from (^9) 



The attachment point of the equivalent mechanical oscillator giving the 
proper moment is presented in Figure 2J. It is observed from (43) that the 
net moment imposed by liquid sloshing in a spherical tank is zero. Heuristi- 
cally, this is easily confirmed by noting that the moment of pressure forces 
acting anywhere in the tank is always zero because an elemental area at any 
point in the tank is normal to a radius vector from the center of the tank. 

This is not true, however, in spheroidal tanks of nonzero eccentricity. It 
may be observed that a positive pressure acting on an element of the tank 
wall to the right and below the tank's center produces a generally negative 
moment (i.e. one tending to turn the tank in the clockwise direction). 
Negative pressures characteristic of the lower Bond numbers will produce 
generally positive moments, hence the positive value of z^ noted at the 
smaller liquid volumes. Increasing liquid volume allows a greater hydro- 
static pressure which has a greater negative moment. Thus for larger liquid 
volumes, the net moment will be negative and the value of z^ also negative. 
The discontinuity between negative and positive values of z^ corresponds 
to the point where the lateral force in Figure 25 is zero. 
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Two- Circle Meniscus Nomenclature 
(b) 


Dimensionless Tank and Liqixid Geometry 
Figure 1 
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Ftriodlc Puloe Function 


Square Wave and Periodic Pulse Perturbing Accelerations 

Figure 3 
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Elgeimode Shapes 
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.genmode Shapes for = 100, V 
Figure 17 
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Shapes for V = 0. 5, e =* 0.5; ®ct “ ^ ~ 

Figure l8 
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Lateral Force Action Point for First Mode Sloshing in Spheroidal Tanks 
with Eccentricities 0.5^ 0.68, and 0.8 
Figure 27 
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APPENDIX A 
SIMBOL LIST 


Quantities are nondimensional unless otherwise designated. When appropriate, 
the relation between a dimensionless variable and the physical dimensional 
one, which is topped by a bar, is given. Underlined variables are generally 
vectors. 


English Alphabet 
a 




(m) 


’tr 

tr 


b 
B 
B. 

A 

B. 
B 

a 

C 

Cl 

''m 

C 

u 


B, 


-Dimensional semi-major axis of container. The character- 
istic length with respect to which the variables are made 
nondimensional. 

-Diagonal matrix whose entries are areas of zones on free 
surface ; 

-Matrix in generalized eigenvalue- vector problem. 

-Coefficient of in Fourier expansion of perturbed 

velocity potential. 

-Coefficient of cos(m(»^t) in Fourier expansion of 
perturbed velocity potential. 

-Semi-minor container axis = b/a. 

-Matrix in generalized eigenvalue- vector problem. 

O 

-Transverse time- varying Bond number = Pg. a /a. 

-Amplitude of B. . 

TfV 2 

-Axial Bond number = pg^a /a (assumed non-negative) . 

-Curvature of meridian of equilibrium free-surface. 

-Contribution of lower intersection to F . 

X 

-Coefficient of sin(nrD^t) in Fourier expansion of periodic 
perturbation. 

-Contribution of upper intersection to F . 

X 

-Matrix approximating the Laplace equation within the liquid 
with zero normal derivative on the equilibrium free-surface 
and w and zero potential on the center line. 

-Foxxrier coefficient in the expansion of R in terms of 
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e 



G 

J4 

h. 

1 

H 

k 


K 

k 

L 

m 




ll4 


-Eccentricity of an ellipse. 

-Unit vector in T1 dii^ction. 

-Unit vector in r direction. 

-Unit vector in s direction. 

-Unit vector in 9 direction. 

-Unit vector in t direction. 

-Arbitrary function of time = L<7/(pa)] f. 

-A generic function of convenience. 

-Lateral force = F^/aa. 

-i'th component of G in surface polar normal coordinates. 
-Dimensional time-varying lateral acceleration (positive 
when acting in direction of increasing x). 

-Dimensional steady axial acceleration (positive when acting 
downward) . 

-A vector function. 

-Mean curvature of free surface = a.J4, 

-Signed scalar magnitude of derivative of position vector 
with respect to i"^^ coordinate direction, i = 1,2,3* 

-See capital eta. 

-Modulus (k = m, where m is the parameter of a Jacobian 
elliptic function). 

-Real quarter period of a Jacobian elliptic function. 

-Unit vector in z direction. 

-Left triangular factor of the matrix S. 

-Index, particularly in Fourier sine expansion of perturbing 
acceleration. 

-Moment of lateral force = M /oa . 

T 

-Equivalent mass for mechanical analogue for kth mode = 

( l+%)^j5.(a2p ) . 

-Unit vector normal to free surface. 

-Unit vector normal to container wall. 

-Pressure = pa/a. 

-Gas pressure = p^a/a. 
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r 


s 

s 

max 

S 

sn 

t 



T 


U 

V 

V 


-Static liquid pressure at a fixed point on equilibrium 
free surface = p^a/a. 

Distance from (R,Z) to upper branch of ellipse along 
r = R. 

-P = max(P:0 ^s^s-6, 0<6< As). 

-Distance from (R,Z) to lower branch of ellipse along 
r = R; 

-Coefficient of in the free- surface boundary condition; 

-A function of convenience. 

-Radial coordinate = r/a. 

-Radial coordinate of point on equilibrium free- surface 
meridian = R/a; 

-Right triangular factor of the matrix S. 

-Principal radii of curvature of equilibrium free- surface 
= R^/a, Rg/a. 

-Radius vector farom origin of tank fixed coordinate system 
= r/a. 

-Arc length along equilibrium free surface meridian = s/a. 

-Symbol for s on plotted output. 

1 - 1-1 

-Tridiagonal matrix S = ^ A TA . 

-A Jacobian elliptic function. 

-Time = [(l+B^)a/(pa^)]^/^ t; 

-Dummy variable. 

-s- coordinates of mid- points of two adjacent mesh intervals. 

-Tridiagonal matrix approximating the free- surface boundary 
operator L; 

-As superscript denotes transposition; 

-Period of = [(l+B^)CT/(pa2)]^/^ T. 

-R^ in system of first order differential equations for 
meridian of equilibrium free- surface. 

-Fluid velocity = Vf = [ (l+B^)cr/(pa)]"^/^ v. 

-Volume = V/(4rTa^/3). 
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Vk 

W 


w 


X 

A 

X 

X 


1 



Y 


z 

z 

zU) 

Z 

Greek Alphabet 
6 

Sjk 

As 

At 

At 

C 

H 

e 

0 


ll6 


-Energy amplitude of k"*^^ mode = Vj^Coa^). 

-Denotes container wall; 

-Complex variable. 

-Weight in integration formula. 

-Z^ in system of first order differential eQ.uations for 
meridian of tank wall* 

-Cartesian coordinate = r cos9 = x/a 

-Maximum lateral mass displacement of mechanical analogue 
= x/a. 

-Radial coordinate of point on container wall meridian 
= x/a (measured from polar axis). 

-Vector of approximations to $ for mesh points not on 
free- surface. 

-Improved approximation for 

-Axial coordinate of point on container wall meridian 
= y/a (measured from equatorial plane). 

-Axial coordinate = z/a. 

-Vector of approximations to § on the free- surface. 

-Axial coordinate of action point of mechanical analog 
= Zj^/a (measured from the tank center). 

-Initial guess for z. 

-Improved approximation for _z. 

-Axial coordinate of point on equilibrium free surface 
meridian = 7^ jo,. 

-Ratio of periodic pulse width to half period to the per- 
turbing acceleration. 

-Kronecker delta =1 if j = k, otherwise = 0. 

-Mesh spacing on s . 

-Perturbing pulse width. 

-Mesh spacing on w. 

-A complex variable. 

-Normal coordinate = T)/a. 

-Normal coordinate of free surface = H/a. 

-Angular coordinate 

-Contact angle. 
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-Spring constant for mechanical analogue for mode 

= 

-Pressure difference across a fixed point on the equili- 
brium free surface = p -p ; 

g o’ 2 

-Approximation to the eigenvalue cu • 

-Guessed \ in Wielandt inverse- iteration. 

-Improved X in Wielandt Inverse- Iteration. 

-Integer index. 

-Integer index. 

-Dimensional liquid density. 

-Dimensional surface tension. 

-Arc length along container wall meridian = x/a. 

-The length of a container meridian from pole to pole 
= Tj-Za. 

-Approximation to I at the jth mesh point. 

-Vector of approximation to I. 

-Vector component of ^ points on equilibrium free- 

surf ace . 

-Vector component of 9 . for points not on equilibrium 
free- surface. 

-Velocity potential = [(l+B^)aa/p]"^^^ 1. 

-$ for k^^ normal mode. 

-3$/3n. 

-Angle between radial direction and container wall meridian. 
-Angle between radial direction and equilibrium free surface 
meridian. 

-Frequency = [(l+B^)a/(pa^)]"^'^^ uj. 

-U) of k"*^^ normal mode. 

-03 of periodic 

-Angular rotation rate of tank- fixed coordinate system 
= [(l+B^)cr/(pa^)]"^/^ n. 

-Euclidean vector norm 


117 


LOCKHEED MISSILES 6c SPACE COMPANY 


Combining Subscripts 

J -Index, usually for points and values belonging to boundaries. 

k -Index denoting noimal mode. 

Sj -Denotes lower circle of intersection, if any, of eq[uilibrium 

free- surface and container wall. 

m -Index, particularly in Fourier sine expansion of perturbing 

acceleration; 

-Also denotes mean or maximum. 

n -Denotes differentiation in the exterior normal direction. 

tr -Denotes transverse (lateral) direction. 

u -Denotes uppermost intersection of equilibrium free- surface 

and container wall. 

w -Denotes container wall. 

a -Denotes the axial direction. 
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